# ------------------------------------------------------------------------------
# Predicts Yield, Test using demographic variables
# Author: Cassidy Shubatt <cshubatt@gmail.com>
# To run: bash 03_demographic_risk_tbl.sh
# ------------------------------------------------------------------------------

# Libraries --------------------------------------------------------------------
message("Loading libraries...")
library(here)
library(yaml)
library(tidyverse)
library(glue)
library(testit) # assert()
library(lfe) # felm()
library(broom) # tidy()
library(data.table) #setnames()
library(stargazer)

temp <- here("code", "07_other_physician_errors", "temp")
u <- modules::use(here("lib", "util.R"))

# Load Data --------------------------------------------------------------------
message("Loading data...")
overnight_lab <- ""
paths <- read_yaml(here::here("lib", "filepaths.yml"))
dems <- readRDS(paths$analysis$demographics)
cohort <- readRDS(glue(paths$analysis$test_cohort)) %>%
  filter(!exclude) %>%
  u$safe_left_join(dems) %>%
  mutate(age = factor(age)) %>%
  setnames(
    c("p__ensemble__stent_or_cabg_010_day__tested", "tile_stent_or_cabg_010_tested"),
    c("yhat", "yhat_decile")
  ) %>%
  mutate(yhat_decile = factor(yhat_decile) %>% relevel(ref = 10))

# Prep regressions -------------------------------------------------------------
message("Preparing demographic regressions...")
dem_vars <- c(
  "age_at_admit", "race_black", "race_hispanic", "race_other",
  "sex_female", "agi_under_25K"
)
covariates <- c("yhat", dem_vars, "1 | 0 | 0 | ptid")
dem_regs <- list()
outcome_vars <- c("test_010_day", "stent_or_cabg_010_day")#, "macetrop_030_pos")

for(outcome in outcome_vars){
  message("Running regression for ", outcome, "...")
  reg_df <- copy(cohort)
  if(outcome == "stent_or_cabg_010_day"){
    reg_df <- filter(reg_df, test_010_day)
  }else if(grepl("mace", outcome)){
    reg_df <- filter(reg_df, !test_010_day)
  }

  form <- reformulate(
    response = outcome, termlabels = covariates
  )

  fit <- felm(form, reg_df)
  dem_regs <- c(dem_regs, list(fit))
}

# Save stargazer ---------------------------------------------------------------
message("Writing regression table with stargazer...")
clean_tbl <- stargazer(
  dem_regs, omit.stat = c("f", "ser", "ll", "aic"), no.space = TRUE,
  dep.var.caption = ""
)

save_fp <- file.path(temp, "demographic_reg_tbl.tex")
write(clean_tbl, save_fp, append = FALSE)

message("Done.")
